Combined Sampled Parameters

Let’s combine and organize our sampled parameters in a data frame. We’ll use an equation-specific function to ensure completeness and repeatability for multi-well matching. This function will apply different methods based on the class of the input parameters list.

# The sampleModel function will sample and combine all parameters with different methods
# based on the parameter set class: MH, MS, etc.
# First define the generic function. Generic function simply calls method:
sampleModel <- function(inputs = c("list of model parameter inputs"),
                        iterations = c("no. iterations for sampling")){
  UseMethod("sampleModel")}

# Next define the method for a MH model
sampleModel.MH <- function(inputs,n){
  # For each MH input parameter, call sample.any function for defineDistribution class list
  qi <- sample.any(n,inputs$qi)
  b <- sample.any(n,inputs$b)
  di <- sample.any(n,inputs$di)
  dlim <- sample.any(n,inputs$dlim)
  
  # Combine into dataframe and return list with class assigned
  return(structure(list(data.frame(qi,b,di,dlim)),class="MH"))
}

# The method for an MS model will be defined in the background for use later

# Call the sampleModel function and MH method to populate a dataframe with all samples
MH.sampled.parameters <- sampleModel(oil.inputs,10000)

# Return in table, displaying first 10 rows only
kable(MH.sampled.parameters[[1]][c(1:10),], align = 'c') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive",
                                      full_width = F, position = "center"))

Historical Rate Synthesis

Using the combined parameter dataframe and the modified hyperbolic equation, we can define a function to use each row of the parameter dataframe as a Monte Carlo iteration and substitute it into the arps equation across the time interval for the sample well. The function will create and return a mcIteration object for each iteration. This object is a list with slots for Iteration name, Parameters, and Synth. Some additional slots will be created and populated later as we approach cost and forecasting.

# A generic function synthesize will be used to apply the correct method based on the 
# class of model used. The class of this model will be in the input list containing 
# the dataframe
# First define the generic:
synthesize <- function(parameters,q.abd,t){
  UseMethod("synthesize")
}

# Now define the MH method of synthesize to accept parameters of class MH, and return a list
# of synthesized rates, and input parameters
synthesize.MH <- function(parameters,q.abd,t){
  # Define builder function to call in apply to build mciteration class from dataframe
  mciteration <- function(x,q.abd,t){
    # Parse dataframe row and pass to MH
    fc <- MH(x[["qi"]],x[["b"]],x[["di"]],x[["dlim"]],q.abd,t)
    # Simplfy parms into string
    parm_string <- paste(names(x),round(x,digits = 2),sep="=",collapse = " ")
    # Combine into mciteration object and return
    return(list("Iteration" = parm_string,"Parameters" = x, "Synth" = fc))
  }
  # Call mciteration in apply across the rows of parameters dataframe and return the result
  return(apply(parameters[[1]],1,mciteration,q.abd = q.abd,t = t))
}

To properly synthesize production data for the sample well, we need to create a time vector from the peak oil rate to the last producing month.

# Define a time interval of actuals to synthesize across, using 100/05-17-081-16W6/00
uwi<-"100/05-17-081-16W6/00"

# Subset production including only 100/05-17-081-16W6/00
prod <- formated_production$CD[which(formated_production$CD$Unique.Well.ID == 
                                              uwi),]

# Find peak oil month for this uwi
uwi.pm <- matrix[["Oil.Peak.Month"]][which(matrix$eachUWI==uwi)]

# Refine production datatable to only include data after peak oil rate
prod_ap<-prod[-c(1:(uwi.pm-1)),]

# Find months on prod
mop <- sorted_matrix[which(sorted_matrix$eachUWI == uwi),"Months.on.Prod"]

# Create time vector starting at 0 for peak oil month and using cumulative sum for days
synth.time <- cumsum(c(0,head(prod_ap$Days.In.Month,-1)))

# Call the synthesize function on the sampled parameters to apply the MH method
# It will return a list of rates, one slot for each parameter sample set
oil.synth <- synthesize(MH.sampled.parameters,q.abd = o.abd, t = synth.time)

# Print first object of oil.synth to show mcIteration object
print(oil.synth[[1]])

## $Iteration
## [1] "qi=460.36 b=0.54 di=0.89 dlim=0.07"
## 
## $Parameters
##          qi           b          di        dlim 
## 460.3625692   0.5439358   0.8949809   0.0700000 
## 
## $Synth
##  [1] 460.36257 330.41732 247.35453 192.62467 155.60928 127.75619 107.50162
##  [8]  91.34834  78.65800  69.38956  60.96063  54.21876  48.38986  43.62112
## [15]  39.41518  35.80406  32.77353  30.03893  27.71393  25.59101  23.70910
## [22]  22.13517  20.62432  19.30858  18.08037

Each mcIteration list, defining a Monte Carlo iteration and historical synthesis, can now be evaluated against the actual data to search for an appropriate model.

Soft L1

Soft L1 loss is a robust, continuous variation of the common Huber Loss that gives less weight to outliers than a L2 loss, making it useful for fitting time series production data. Our Soft L1 definition will return the individual loss of each point and an average cost for all residuals input.

soft_L1<-function(r){
  IndvCost <- 2*((1 + abs(r))^0.5 - 1)
  Cost <- sum(IndvCost,na.rm = TRUE) / length(IndvCost)
  return(list(IndvCost,Cost))
}

To apply the loss equation we’ll define a function to apply the Soft L1 loss to the residuals of each model match. Due to the non-linear nature of production data, the loss function will be applied to the log residuals. The function will populate the IndvCost and Cost slots of each mcIteration object. For this illustration, the function will be defined in the background and we can see the result for the first iteration object.

## $Iteration
## [1] "qi=460.36 b=0.54 di=0.89 dlim=0.07"
## 
## $Parameters
##          qi           b          di        dlim 
## 460.3625692   0.5439358   0.8949809   0.0700000 
## 
## $Synth
##  [1] 460.36257 330.41732 247.35453 192.62467 155.60928 127.75619 107.50162
##  [8]  91.34834  78.65800  69.38956  60.96063  54.21876  48.38986  43.62112
## [15]  39.41518  35.80406  32.77353  30.03893  27.71393  25.59101  23.70910
## [22]  22.13517  20.62432  19.30858  18.08037
## 
## $IndvCost
##  [1] 0.3134137610 0.2965712496 0.3087135803 0.3350164450 0.2721060015
##  [6] 0.2146343907 0.0235178085 0.0433602151 0.0245615456 0.0119766911
## [11] 0.0008584473 0.0732817531 0.0628184311 0.0506526532 0.1212306678
## [16] 0.1705506027 0.1206529242 0.1697310024 0.2106775721 0.2816378928
## [21] 0.2649857987 0.2121464179 0.2548659827 0.2265560757 0.2848290830
## 
## $Cost
## [1] 0.1739739

Using the cost function applied to each iteration, we can minimize the cost to find our bestfit model. Using this bestfit, we can also visualize the cost model results for each production point in the fit.

Using the log residuals and the Soft L1 cost function, let’s find the top 100 models from our initial 10,000 iterations and see how the matches vary from the bestfit. We’ll colour them by average cost to compare the forecasts.

Top 1% of Models Using Soft L1 Loss

We can see that the top 100 forecasts do a good job of covering the range of actual oil rates, with the lowest cost forecasts tightly grouped around the actuals. However, one observation is that the forecasts are not evenly distributed above and below the actuals to appropriately capture uncertainty in the forecasts. Let’s explore another loss function to see if we can improve this behaviour.

PDCA Loss

The Soft L1 loss function does a decent job at finding a best fit model, but is there a better loss function that can be used? Let’s look at a loss function defined in SPE-174784-PA (ADD LINK HERE AS REFERENCE) by Fulford et al. that takes a similar approach, but is designed to capture more uncertainty in the production data. This loss function uses the average and standard deviation of the individual cost values across the fit, comparing them to a baseline fit. This will prove particularily useful when we move into the PDCA methods later in the CAST program.

For now, let’s define this loss function as PDCA Loss and see how it improves the fits from our Monte Carlo.

# As our Soft L1 was defined, pdca_loss will return slots for IndvCost and Cost
pdca_loss<-function(r, Eps_mean_min = 0, Eps_std_min = 0){
  IndvCost <- r^2 # pdca loss simply uses L2 loss for individual cost
  
  Eps_mean <- mean(r,na.rm = TRUE) #Mean of residuals
  Eps_std <- sd(r,na.rm=TRUE) #Std of residuals
  
  Cost <- ((Eps_mean - Eps_mean_min)/.1)^2+((Eps_std - Eps_std_min)/.01)^2
  return(list(IndvCost,Cost))
}

Just as we did with the Soft L1, let’s use the PDCA loss function to find the top 100 forecasts from the same 10,000 iterations. We’ll set the minimum mean and minimum standard deviation of epsilon to 0, assuming a perfect match. We’ll colour by average cost and compare to the Soft L1 method.

Comparing Top 1% of Models Using Soft L1 and PDCA Loss

We can see that the PDCA Loss function has improved the distribution of fits above and below the actuals, a behaviour that we would like to carry forward because it helps to describe a high and low fit, as well as a bestfit. Lastly, let’s look at the parameter distributions for the top 1% of models to get an idea of the solution space using the PDCA Loss method.

Parameters of P1 model matches using PDCA Loss method

Summarize observations from these distributions.. They are mostly lognormal. Therefore, do we need to step into MCMC to get a proper PDCA distribution going forward? Keep these in mind and compare later on..

Summarize findings in these two charts and how we’re going to carry PDCA loss going forward in the rest of the evaluation / process as it provides a better probabilistic representation. Show distribution fits of the accepted parameters? Or do this at a later time? Or here just FYI??

Burn-In

The burn-in process is an initial convergence step that allows the MCMC to start from an equilibirum position, in this case the bestfit from the 10,000 Monte Carlo runs. While the burn-in process is often described as not needed or biased, in the case of PDCA we would like maintain its functionality to give us a baseline fit for evaluation of the PDCA loss function. In plain words, there is more uncertainty in production data that is hard to fit a curve to rather than data that is easy to fit. In this way, we can honour that higher level of uncertainty when needed.

Metropolis-Hastings Algorithm

With the bestfit defined as the starting point, we can start sampling through the solution space to carry out the MCMC process. The Metropolis-Hastings algorithm is a method for obtaining a sequence of random samples from multi-dimensional distributions like we have in PDCA. Simply defined, the Metropolis-Hastings method starts at an arbitrary multi-dimensional point (the start of the Markov Chain) and randomly chooses the next multi-dimensional point to evaluate given jumping distributions for each dimension. In our case the inputs parameters for the MH ARPs equation make up each dimension. For each iteration, the candidate sample is compared to the previous sample using an acceptance ratio. If this ratio is greater than 1, the jump is automatically accepted as being more probable than the previous sample. If the ratio is less than 1, the point is less probable than the previous sample, but can still be accepted. Using a random a random probability, sometimes the move will be accepted and sometimes it will be rejected, with acceptance ratios closer to 1 being given a higher probability of acceptance. Thus, we will tend to stay in high density regions and define the desired distribution after sufficient iterations.

The jumping distributions we will use are defined as follows:

Use some formulas to describe the N(i-1,SD) of each jumping distribution

Applying these jumping distributions for 10,000 iterations on our sample well, we can follow the chain of accepted or rejected iterations as the posterior distribution is defined.

Accepted input parameters using MCMC method

The resulting accepted distributions are dependent on the jumping distributions, which need to be tuned appropriately, but we can see the posterior distributions take shape. With more samples than the Monte Carlo illustration, each distribution is well defined, following normal or lognormal distributions.

Add note about staying in the proper solution space (prior distributions) through out entire MCMC process - this is crucial to have the equations make physical sense

Results

Let’s compare the top Monte Carlo matches to the MCMC matches to see how the accepted solution space matches the oil rate behaviour.

Comparing Top 1% Monte Carlo Matches to MCMC Matches

We can see that the matches have a tighter spread around the actuals, returning more forecasts within the acceptable solution space. The MCMC and Metropolis-Hastings algorithm method has achieved the desired result, returning thousands of acceptable models that define the uncertainty in forecasting the sample well. Using the same methods derived earlier in the markdown, let’s compile the accepted forecasts into P10, Mean, P50, and P90 durves to summarize the range of outcomes.

Finally, let’s summarize the EUR distribution and the EURs represented by the probabilistic forecasts.

PDCA Results Comparison

MCMC results comparison of MH (left) and MS (right) models

MCMC EUR comparison of MH (left) and MS (right) models